Influence of conservative corrections on parameter estimation for EMRIs 



OO 

o 



X 



PACS numbers: 



E. A. Huerta and Jonathan R. Gair 
Institute of Astronomy, Madingley Road, CB3 OHA Cambridge, UI^ 
(Dated: April 22, 2009) 

We present an improved numerical kludge waveform model for circular, equatorial extreme-mass- 
ratio inspirals (EMRIs). The model is based on true Kerr geodesies, augmented by radiative self- 
force corrections derived from perturbative calculations, and in this paper for the first time we 
include conservative self-force corrections that we derive by comparison to post-Newtonian results. 
We present results of a Monte Carlo simulation of parameter estimation errors computed using the 
Fisher Matrix and also assess the theoretical errors that would arise form omitting the conservative 
correction terms we include here. We present results for three different types of system, namely the 
inspirals of black holes, neutron stars or white dwarfs into a supermassive black hole (SMBH). The 
Q.^ \ analysis shows that for a typical source (a IQMq compact object captured by a 10^ Mq SMBH at 

, a signal to noise ratio of 30) we expect to determine the two masses to within a fractional error 

' of ~ 10"*, measure the spin parameter q to ^ 10"*"'' and determine the location of the source on 

I the sky and the spin orientation to within 10""^ steradians. We show that, for this kludge model, 

, ^ , omitting the conservative corrections leads to a small error over much of the parameter space, i.e., 

^ • the ratio TZ of the theoretical model error to the Fisher Matrix error is 7?. < 1 for all ten parameters 

' in the model. For the few systems with larger errors typically TZ < ?> and hence the conservative 

, corrections can be marginally ignored. In addition, we use our model and first order self-force 

fSJ ■ results for Schwarzschild black holes to estimate the error that arises from omitting the second- 

^Nj ' order radiative piece of the self-force. This indicates that it may not be necessary to go beyond first 

, order to recover accurate parameter estimates. 

cr. 

^ ; I. INTRODUCTION 

> : 

OO , One of the most important sources of gravitational waves (GWs) for the Laser Interferometer Space Antenna 
(LISA) will be the extreme-mass-ratio inspirals (EMRIs) of stellar mass compact objects into massive black holes 
^ (BHs). These systems occur in the centres of galaxies where massive black holes are surrounded by clusters of stars 
■ which contain large numbers of compact stellar remnants. The heavier black holes sink toward the centre as a result 
' of mass segregation during random encounters among stars. These encounters will gradually perturb the orbits until 
one of the stellar mass black holes passes close to the central black hole. If the impact parameter is sufficiently small, 
enough orbital energy will be radiated in gravitational waves to leave the remnant bound to the BH. Thereafter, the 
small black hole orbit gradually decays via GW emission. Eventually the black hole may end up on an orbit where it 
is radiating GWs continuously in the LISA band. Even then, the black hole may complete several hundred thousand 
orbits before crossing the horizon. 

The GWs emitted by these sources will be rich with information. They will provide a map of the space-time exterior 
d ' to the large body and the response of the horizon to tidal forces. Furthermore, tests of the no-hair theorem will be 
possible as LISA will be able to measure the black hole's mass, spin and quadrupole moment to fractional accuracies of 
10"'^ [l[. Testing this theorem will allow us for the first time to confirm or falsify the predictions that general relativity 
makes about black hole solutions. On the other hand, we could discover non-black hole systems in nature, e.g., boson 
stars or naked singularities. Additionally, the accurate system masses and spins that LISA should determine will be 
of great importance for understanding the stellar populations in the central parsecs of galactic nuclei. 

In contrast to comparable mass binary systems, where modelling may be done using post-Newtonian (PN) theory 
in the early inspiral and numerical relativity for the final orbits and merger, EMRI GWs may be obtained using 
Black Hole Perturbation Theory (BHPT). This is because the extreme mass ratio, r] = m/M, can be used as a small 
expansion parameter. Such EMRIs will be visible to LISA for central black holes with masses in the range 10'* — 10^ 
Mq and out to redshifts z « 1 ^, i3||. As we will show later, the small mass ratio rj has an important effect on the 
number of orbits that the small body performs before merging. 
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Detection of EMRIs will be difficult, since it is expected that the signal from EMRIs will be about an order of 
magnitude below LISA's projected instrumental noise. In fact, at lower frequencies, below 2 mHz, the signal will be 
orders of magnitude below the confusion noise due to unresolved galactic binaries. However, signals can be identified 
in noisy data using matched filtering. This requires the use of theoretical templates whose phase remains accurate to 
one cycle over the r]~^ cycles of the waveform generated while the orbit of the stellar mass compact object is in the 
strong curvature region, close to the large black hole. Hence, there is a pressing need for accurate theoretical models 
of the waveforms. Modelling this emission is a difficult problem, although there has been great progress in recent 
years. 

BHPT relies upon the fact that the stellar mass compact object can be treated as a small perturbation of the 
gravitational field of the central black hole. Over short timescales, the compact object follows a Kerr space-time 
geodesic since we can neglect back-reaction for that short stretch of time. A Kerr geodesic is characterized by three 
constants of motion, namely the energy, E, the angular momentum about the hole's spin axis L^, and the Carter 
constant Q, which is a relativistic generalization of the third integral of motion, which arises from the separation of 
the equations of motion for orbits in an axisymmetric gravitational potential. In the spherical limit, i.e., for a non- 
rotating central black hole, Q reduces to the square of the angular momentum projected into the equatorial plane. 
Over longer timescales, ~ M /"q, radiation-reaction causes the orbit to evolve adiabatically. This can be characterized 
by changes in the geodesic orbital elements and Q. As the compact object inspirals, it eventually passes through 

an innermost stable orbit where adiabaticity breaks down. It then follows a geodesic plunge orbit and is swallowed 
by the central black hole. 

First order radiative BHPT in a Kerr background is described by the Teukolsky formalism. Unfortunately, 
Teukolsky-based waveforms are very computationally expensive to generate. However, there exist various fami- 
lies of approximate waveforms that capture the main features of true signals and are much quicker to generate. One 
family are the post-Newtonian (FN) waveforms, which are both analytic and easy to generate. Recently, a class of 
FN waveforms were developed for the EMRI case by Barack and Cutler ^] . These GWs were constructed considering 
the lowest-order quadrupole waveforms for eccentric binaries on Keplerian orbits derived by Peters and Matthews [H] , 
but the orbits were corrected to include the effects of pericenter precession, Lense-Thirring precession, and inspiral 
due to radiation reaction. This family of GWs are useful for scoping out data analysis issues, e.g., computing the 
Fisher information matrix to estimate parameter measurement accuracies. However, a small black hole in a close 
orbit about a much larger one is too relativistic a system for FN analyses to be valid. The FN expansion is unlikely 
to be reliable in the EMRI limit since the object spends so much time in the strong field region of the space-time 
where v ^ c. But, it is in this zone where most of the GWs observable to LISA are generated. In order to bridge this 
gap, we again rely upon BHFT. 

Another family of EMRIs waveforms are referred to as "numerical kludge" (NK)_6j. These are constructed by 
combining an exact particle trajectory with a fiat space-time wave-emission formula. This framework also aims to 
capture the main features of the waveform accurately. Comparison to Teukolsky-based waveforms indicates that NK 
waveforms for geodesic orbits are quite faithful as the overlap between the two types of waveforms for geodesic orbits is 
greater than 0.95 over a considerable portion of the inspiral parameter space @. However these models are incomplete 
in their treatment of the orbital evolution under the self-force since the model currently does not include conservative 
radiation-reaction terms, and the waveform-emission formula and phase space trajectories are approximate. 

There is currently some debate in the literature 0, [1] as to the importance of the conservative pieces of the self- 
force for waveform modelling, both in terms of source detection, and parameter estimation. The purpose of this paper 
is to shed some light on this issue using the numerical kludge waveform family. We will describe how to include 
conservative corrections in our model, and then estimate their infiuence on parameter extraction using the formalism 
recently developed by Cutler and Vallisneri 

This paper is organized as follows. In section[TT]we construct asymptotic observables for our kludge waveform model, 
namely the orbital frequency and its first time derivative. We show how we can effectively include both components 
of the self-force: radiative and conservative, to second and first order, respectively, by comparison to FN waveforms. 
We also discuss the relative importance of the various conservative corrections to the phase evolution of our waveform 
signal. In section Hill we use our waveform model to assess the importance of the second order radiative piece that will 
be missing once first order accurate self-force waveforms are generated. We do this by comparing to self-force data 
obtained by Barack and Sago for circular equatorial orbits around a Schwarzschild black hole. In section HVl we will 
use our inspiral kludge waveforms including modulations from LISA's response function to estimate the noise-induced 
parameter errors using the Fisher Matrix formalism. These provide confirmation and extension of the results obtained 
using the analytic kludge model of Barack and Cutler [4|. In section |V] we employ the recently developed formalism of 
Cutler and Vallisneri [9| to estimate the theoretical parameter errors that would be introduced by waveform template 
inaccuracies. Specifically, we study the importance of the conservative corrections by computing the theoretical errors 
in the parameters that would arise from omitting conservative corrections from the model. Finally, we summarize our 
work in section IVTl 
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II. KLUDGE WAVEFORM 



Our waveform is based on the numerical kludge approach described elsewhere [y, but we summarise the 
formalism here. The aim is to capture the main features of the waveform inspiral in the strong field regime. Because 
an accurate description of the orbital and phase evolution of the small body is not yet available, we try to circumvent 
this problem in various steps. The system's extreme mass ratio m/M ^ 10^^ guarantees that gravitational back- 
reaction effects occur on timescales much longer than any orbital timescale. Hence, we can assume that the inspiralling 
object instantaneously follows a Kerr space-time geodesic. We slowly evolve the parameters of the geodesic using a 
prescription for the rate of change of energy, angular momentum and Carter constant and calculate the resulting 
trajectory of the inspiralling body in the Boyer-Lindquist coordinates of the Kerr space-time of the central black hole. 
The gravitational waveform is then constructed by identifying these coordinates with spherical polar coordinates in a 
pseudo-flat space and using weak-field emission formulae. 

In this paper, we restrict our attention to circular equatorial orbits. For such orbits, the inclination angle and 
eccentricity remain constant, i.e., circular-equatorial orbits remain circular-equatorial [l^ . The energy E and angular 
momentum of such orbits is given in terms of the Boyer-Lindquist radius of the orbit, p, as [Tsj 

E _ l^2{M/p)±{a/M){M/pf^^ 



m 



^Jl-3iM/p)±2ia/M) (M/pf^' 



L. ^ ^ / p y /2 1 T 2 (g/M) {M/pf^ + {a/uf {M/pY 

^1-3 (M/p) ± 2 (o/M) {M/pf^ 

where the upper (lower) sign is for prograde (retrograde) orbits and a stands for the spin of the central black hole. 
From now on we will write q — a/M . The Carter constant, (5 = 0, for equatorial orbits. 

To obtain the orbital evolution of the compact object, we need to calculate the evolution of the energy E and angular 
momentum L^. We can evaluate these quantities by equating their rate of change with the fiux carried away by the 
gravitational waves, E and L^- We will use the radiation fluxes derived by Gair & Glampedakis [ll| to compute the 
compact object's orbital evolution. These fluxes are based on 2.5PN expressions derived by Tagoshi [ij], augmented 
by fits to accurate BHPT results and with various consistency conditions imposed to ensure physical behaviour for 
near-circular and near-polar orbits. 

To evolve a circular orbit we need only to specify the angular momentum or the energy fiux, as they are related by 
the 'circular goes to circular' rule 

E{p) = ± Up) = f7(p)i.(p), (2) 

where d0/dt = f7(p), is the azimuthal velocity of the orbit. We choose to work with the evolution of L^, which has 
the following form [lli] 

44711 /A// \^ 33 ^[mV , , , , n ] 

~ "9072" Vj) ^Te'^iyj ^ ^ ^ Teukolsky fits ^ . (3) 

The 'higher order Teukolsky fits' are given in [ll|. We do not give these explicitly here, as they are not needed to 
derive the conservative corrections. However, we will include them to evolve orbits when we generate waveforms. 
The evolution in time of the radial coordinate is given by 

Using the exact geodesic expression for dp/dLz generates inspirals that are closer to Teukolsky based evolutions than 
expanding the above expression at 2PN order. However, we will need the 2PN expression in the following, which is 

dp 64 f mV \ ^ 743 fM\ / 133 \ /M^^/^ 



mM 



/34103 , 81 2\ ( ^ " 



V 18144 16 



• (5) 
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As before i] = m/M . Up to this point our analysis has been incomplete as we have only considered the three integrals 
of motion. There are also three positional constants of the motion, which basically label the position of the test particle 
along the geodesic trajectory at some fiducial time. The evolution of these is non-trivial because the self-interaction 
between the inspiralling body and the central black has two main pieces. The only one we have considered so far is 
the dissipative or radiative self-force and affects the evolution of the orbit. However, in real inspirals there is a second 
piece called the conservative self-force, which affects the principal constants of the motion. In other words, it changes 
the frequency of an orbit at a given radius but does not cause the orbit to evolve. We need to include this piece of 
the self-force as it could lead to several cycles of phase discrepancy in our kludge waveform over the inspiral. 

The conservative self-force has two parts. One part is oscillatory and averages to zero, whereas the other piece 
accumulates and affects the phasing of the waveform over time. In our kludge we may account for this effect by 
changing the (j) frequency. Specifically, we include this effect by writing 



This equation includes the phase derivative for a geodesic, labeled by the subscript "geo" , and a frequency shift 
which will depend on the instantaneous orbital parameters. A problem arises here because to compute the necessary 
frequency shifts within our framework, i.e., BHPT, would require self-force calculations. At present, these are known 
only for particles moving on circular geodesic orbits around a Schwarzschild black hole, and only in a particular gauge 
[lO| . But, there is an ongoing effort to extend such calculations to generic inspiral orbits in a Kerr background. The 
main challenge in this effort has to do with gauge freedoms as we will discuss later. 

However, we do know conservative corrections in the post-Newtonian framework up to 2PN order which include 
spin-orbit, spin-spin effects and finite mass contributions [15]. We can combine these expressions with the radiative 
self-force obtained by Tanaka et al. |16|| based on the Teukolsky and Sasaki-Nakamura formalisms for perturbations 
around a Kerr BIT which includes terms of order . Equipped with these results we will extend the method proposed 
by Babak et al. [y| to include conservative corrections in the kludge model. They computed the IPN conservative 
correction for circular orbits in the Schwarzschild space-time, and we now extend that calculation to derive 2PN 
conservative corrections for circular orbits in the Kerr space-time. 

The idea is to correct the kludge expressed in a particular coordinate system, in order to ensure that asymptotic 
observables are consistent with post-Newtonian results in the weak field. In particular, we aim to modify the orbital 
frequency and its first time derivative. By modifying these two quantities, we can both identify coordinates between 
the two formalisms and find the missing conservative pieces. 

From equation ([J), at 2PN order the orbital frequency takes the form 

which we now augment by including conservative corrections 
A4> 




+ 77(^^0 + di(^y)+(di.5 + 9/1.5) ( — ) +d2[ — ] ]\. (8) 

We will use this expansion for figoo only to derive the conservative corrections. As with dp/dL^, it has been shown 
that more reliable waveforms can be obtained by including the full geodesic frequency where it is known, c.f., ([5]), and 
this will be the approach used in section [TVl It is inconsistent to include some pieces of the evolution at arbitrary PN 
order, while including the conservative corrections only at 2PN. However, the lower order effects that we are including 
with greater accuracy do have a more significant impact on the waveform. Moreover, it has been found in the past 
that including higher order terms where they are known is the right strategy to obtain accurate waveforms [Tll | . 

It is worth noting that the expansion in p is an expansion in = M/p. To derive the conservative corrections, we 
choose to leave the time derivative of the radial coordinate unchanged and given by equation (H])/®. This amounts 
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to a choice of gauge such that the rj^ piece of dp /At vanishes. Differentiation of ([8]) then gives dVt/dt for the kludge, 

df7 96 ri f mV^'^ { , Mf 743 /5 , 743, 

dF ^ TM^Uj \' + '''^° + 7r536+H3'^'"^ ° 

/MN^^V. 157 A , „, 157, 
Vp"y ( 47r - — g + 77 ( 47rdo + 2di.5 + 9 ( 2/1.5 - — rfo 

M\V34103 81 , /34103 , 81 ,, 3715, 7 



I 1 I — + ri do H Q^dn di H — do 

p J VI8I44 16 'V18144 16 1008 3 



(9) 



We may now write down a coordinate transformation to relate our coordinates with those used in the post-Newtonian 
formalism. 

+ ^(^o+(f)cr+(f)'''(c.5 + ,,..5)+(f)'c.)}, (10) 

where R denotes the post-Newtonian semi-major axis. We can now substitute this expression for the coordinate 
transformation into relations ([8]) and ([9]). 

The final stage of the computation is to compare the expressions for and Q, where a dot denotes d/dt, with the 
available post-Newtonian expansions. The post-Newtonian expansions are available to higher order in the mass ratio 
?7, but we keep rj only to the same order as the kludge, Eq. ([9]). The post-Newtonian expressions for the orbital 
frequency and its first time derivative are given by 

-(^r(-r')}' 

where tox = M + m and L is a unit vector directed along the orbital momentum. Additionally, 
!ip» = f,m»/»„"/=|l-(|2 + ii„)(mT..)V-' + (to-/J)(„,T^) 



/34103 81 o / 13661 



+ 1,18144+16'^ +^ 

where the spin-orbit parameter (3 is given by /? = 1/12 (ll3M^/m^ + 7577) L-q and the constant C will be determined 
from the kludge prescription. Note that we have assumed that the spin of the inspiralling black hole is negligible with 
respect to the central one and so we ignore spin-spin interactions . 

We now rewrite these expressions in a convenient way to take the EMRI limit, by writing ttit — M(l + 77). We find, 

mV /15 169 \ 1 

r) (Y^leV ' ^''^ 
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„„.,,..., , 3 /2591 13571 \ M 
^PN = -— (-) ^^+2-[-336+-67^Vr 

22115 81 , /87044 , /567 A\\/A/^^ 



648 16^ ' V 567 ^ \32 J J \ R 



(14) 



A direct comparison between the two expressions for the orbital frequencies and their first time derivatives aUow us 
to solve simultaneously for bi, 61.5, S1.5, cq, ci, C1.5, gi.5, C2, do, f^i, f^i.s, /1.5, and ^2- We find that the non-vanishing 
parameters are 

1 845 , 1 , 1975 

O1 = 1, Cn = , C-i — , On = — , "1 = 

' "4' 448' 8' 896 

9 27 191 91 

Ci.5--5^, di.5--Y^7r, /i.5--Y^, ^1-5-"^ 

2065193 , 1152343 , , 

Co — , do = . 15) 

677376 ' 451584 ^ ^ 

These parameters not only give us the missing conservative pieces, but also the value of the constant C of equation 
and the coordinate transformation. We find that ( = —243/32 and hence, the first time derivative of the 
post-Newtonian orbital frequency at 2nd order including spin effects and conservative corrections is given by 

■ 96 ?7 3 /2591 13571 \ M 



157 / „ 149 \\ /M\ 
4vr--, + ,fl2.-— 



3/2 



/22115 , 81 2 /87044 81 fM 



■ V 648 i'^^+nw + i^^jjlijl ^''^ 

Having obtained the conservative pieces, we can now assess how important they are in determining the phasing of 
the waveform. First, we will look at the number of cycles that the orbiting stellar-mass compact object performs in 
the last year of inspiral before plunge, and how many of these are contributed by each term in Q. We focus on this 
particular period of time as that is when most of the GWs are radiated in the LISA band. 

In Figure[T]we show the number of cycles in the last year of inspiral including the contributions from IPN, 1.5PN 
and 2PN conservative corrections. We consider six different binary systems which consist of a central Kerr black hole 
of specific mass M = M/Mq = {10^, 10^, 10^} and an orbiting compact object of specific mass m = m/MQ = {1, 10}. 
We consider both prograde and retrograde orbits. 

As mentioned before, even though the GW signals are instantaneously below instrumental and confusion noise, we 
should be able to detect the signals with high SNR observing them for 77"^ cycles, and using matched filtering. But, 
this requires the template waveform to match the true signal to better than half a cycle over the observation. In 
Figure [5] we assess the importance of including the conservative corrections by showing how long all the conservative 
pieces combined take to contribute one cycle to the GW phasing. We consider both prograde and retrograde orbits 
for all of the binary systems shown in Figure [1] 
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FIG. 1: This figure shows the number of gravitational waveform cycles generated by the inspiral of a stellar-mass compact 
object of specific mass ffi falling into a central Kerr black hole of specific mass M over the last year before plunge. The 
number of cycles is derived including conservative corrections up to order 2PN. We show results for both prograde {q > 0) and 
retrograde {q < 0) orbits for the following binary systems, {M, m} = ({10^ 1}, {10^ 10}, {10^ 1}, {10^ 10}, {lO'^, 1}, {10^ 10}), 
from top to bottom, respectively. 
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FIG. 2: This graphic shows the time in years taken for the conservative pieces to contribute one cycle to the GW phasing for a 
given binary system as a function of q. The curve has the same shape for aU binary systems. The vertical scale changes according 
to the binary system under consideration, and is expressed in units of 10" years, where n depends on the system. Hence if 
{M, m, n} defines our system, we show the cases: {10^ 1, 4}, {10^ 10, 3}, {10^ 1, 2}, {10^ 10, 1}, {10^ 1, 0}, {10^ 10, -1}. 

From Figure [2] we learn that for a binary system of 10 -I- 10^ solar masses, the conservative corrections contribute 
one cycle in about 60 days when q = —0.9. In all other cases it takes longer. The next step is to assess the importance 
of the highest order conservative term, i.e., the 2PN term. We do this by playing the same game as above, but 
considering the 2PN term only. This is shown in Figure [31 

10" years 
10 h 



2 ^ 



-0.5 



0.5 



FIG. 3: This graphic shows the time in years taken for the 2PN conservative piece to contribute one cycle to the GW phasing 
for a given binary system as a function of q. The notation is as for Figured and we show the following {M,m, n} cases: 

{lo^ 1, 7}, {lo^ 10, 6}, {lo^ i, 5}, {lo^ 10, 4}, {lo^ 1, 3}, {lo^ 10, 2}. 
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We learn from Figure [3] that the 2PN conservative contribution does not have a significant impact on the number 
of cycles. In fact, referring to the 10 + 10^ binary system considered above, we observe that the 2PN term needs 
about 141 years to contribute one cycle for a value of g = 0.9. In any other case, it takes longer. We also observe 
that in Figure [5] the number of years increases as q increases, whereas Figure [3] shows the opposite behaviour. This 
arises because all of the conservative pieces have a positive contribution except the 1.5PN term. Furthermore, this 
term is the dominant one. Were it not included in the kludge, the number of years taken by the conservative pieces to 
complete a cycle would be much smaller. However, because the 1.5PN term balances and even overcomes the other 
positive contributions we get the behaviour portrayed in Figure^ Roughly speaking, an effect will have to be included 
in a matched filtering template if it contributes of the order of one cycle over the observation. These results suggest 
that it will be necessary to include conservative corrections in our waveform model, but it may not be necessary to 
go beyond 2PN order. 



III. A COMPARISON WITH THE SELF-FORCE FORMALISM FOR NON-SPINNING BLACK HOLES 

The conservative pieces that we have now included in our model, Eq. (jlSp . were obtained in the post-Newtonian 
framework. In the future, we would want to obtain a more accurate solution for the conservative correction by using 
fits to self-force calculations once these are available. However, it is unlikely that results beyond first order in the mass 
ratio will be available before LISA files and we have seen above that second order radiative corrections contribute to 
the phase evolution at the same level as the first order conservative part of the self-force. Our model includes both 
corrections, and so we can use it to assess the relative contribution from the first order conservative self-force and the 
second order dissipative self-force. To do so, we will use the self-force results recently obtained by Barack and Sago 
[l0| for a circular orbit around a Schwarzschild black hole. 

The main complication in self-force calculations comes from gauge freedom — the gauges in which the local and 
global fields are most easily computed are different. Barack et al. [10| found a way to circumvent this problem in the 
non-rotating case by solving both the perturbation equations, and the global retarded field in the Lorenz gauge. In 
our kludge prescription, by choosing only to modify and not dLz/dt when comparing to the PN model, we have 
effectively chosen a gauge in which the ry^ corrections to the radiative part of the self-force vanish. 



A. Assessing self— force calculations using our model 



We can compare self-force results to the kludge results in the same way that we computed the conservative cor- 
rections above, namely by comparing asymptotic observables — the orbital frequency and its first time derivative. 
The second order piece of dil/dt as a function of fl depends on both the first order conservative and the second order 
radiative self-force corrections, and therefore cannot be determined completely from the first order self-force. We 
want to understand whether first order accurate self-force calculations will be sufficient for parameter estimation with 
LISA, by assessing the importance of the second order radiative piece to the orbital phase evolution. 

Following [l^, the orbital frequency of an orbit in the Schwarzschild metric under the influence of the radial 
self-force, Fr, is 



SF 



no 



2Mm 



(17) 



where ilo = \/M/r^, and r stands for the radial coordinate in the self-force framework. 

Barack and Sago [lO] gave a fit to F^, but the terms they proposed were not motivated by PN expansions, 
easier comparison to our results, we can derive an alternative fit to their results of the form 



For 



F^'ir > M) 



M 

flQ + fll 1^ fll 

r 



3/2 



a2 



M 



(18) 



where the various coefficients have the values 



ao = 2.01638, ai = -8.98147, ai.5 = 12.20270, as = -16.51406. (19) 

This fit reproduces the numerical data within the numerical accuracy considered by Barack et al. [lol | in the range 
8M < r < 150M. Note that the leading term F^ ~ 2.016m^/ro ~ 2m? /tq is consistent with the Keplerian SF which 
describes the back-reaction effect that arises from the motion of the black hole around its center of mass. 
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To fit the data in the range 6M < r < 8M we can use the expression 



u M , fM 
bo + bi— + bi.ro — 
r \ r 



3/2 



+ b2 



M 



where the various coefficients are given by 

bo = 2.44194, bi = -29.44529, 61.5 = 89.73089, 62 = -99.32447. 
We can now plug these fits into Eq. p7|) to obtain fisF- 



^SF = ^0 

where the various coefficients have the values 



f , M , fM 
1 + T]\ Ao + Ai— + A1.5 — 
\ r \ r 



3/2 



+ A2 



M 



(20) 



(21) 



(22) 



6M <r<8M : Ao = -1.22097, Ai = 15.94361, Ai.,5 = -44.86544, A2 = 37.38153; (23) 
r>8M : Aq = -1.00819, Ai = 5.49892, A1.5 = -6.10135, A2 = 5.78267. (24) 

To compare the self-force orbital frequency, Eq. (j^ . to the kludge O, Eq. ([5]), we need the coordinate transfor- 
mation that relates the kludge radial coordinate p to that of the self-force formalism, r. We use a transformation of 
the form 



M fM 

p = r< 1 + ax V ai.5 — 

T \ r 



\ r \ r 



3/2 



Oi2 



M 



(25) 



3/2 



We find the non-vanishing coefficients are given by 



6Af < r < 8M 
r > 8M 



j3o = 0.89731, 
/3o = 0.75546, 



/3i 
Pi 



-9.15958, 
-2.19646, 



-/52 



/3l.5 
/3l.5 



M 



24.25543, 
-1.58730, 



P2 - 
P2 



-23.21984 
-2.15393. 



Using the coordinate transformations, Eqs. (|25|) . ([26ll . (pTj) . in Eq. dH) recovers Eq. (fTT]) . 

Under the influence of the conservative self-force, the angular momentum of an orbit undergoes a shift 



L = Lo\l 



2Mm 



(26) 
(27) 



(28) 



where Lq — ^ Mr^ j (r — 3A/). The first order radiative self-force is given by Eq. ([H]), but with p, the kludge 
radial coordinate, replaced by r, the self-force radial coordinate. As before, the evolution of f is best expressed as 
r = LzliArjALz) with 



M- 



dr 

dz: 



2(r/Af - 3)3/2 1 M 



6 + f (.-i-^ + .o + .o.5(^) 



1/2 



M\ 



M\ 



3/2 



where the various coefficients can be found to be 

6A/ < r < 8M : e_i = -2.44194, 
ei = -297.973, 
r > 8M : e_i = -2.01638, 
= -49.5422, 



r ) ^'^^'^ \ r ) 



60 = -14.7936, 
5 = -269.193, 
eo = 3.11681, 
£1 5 = -36.6081, 



£2 



M 



(29) 



eo.5 = 179.462, 
£2 = 595.947; 
eo.5 = 24.4054 
£2 = 99.0843. 



(30) 
(31) 
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By construction, in this gauge, we are missing a contribution from the second order part of L^. However, we have 
already computed a similar expression within the kludge approach, Eq. which includes that contribution. Using 
the coordinate transformation, Eqs ((25l) -([27 l) . we can use the kludge result to find the second order piece of in the 
self-force gauge 

where 6M < r < 8M : Iq 
r>8M: Iq 

The importance of the second order radiative piece of the self-force is readily assessed by computing the difference 
between an evolution using only Eq. ^ for Lz and one that also includes the second order correction above. We 
can now compute the number of cycles that the stellar mass compact object (CO) performs before plunge for various 
binary systems in three separate cases — 1) ignoring all 77^ corrections to f (i.e., setting 6Lz and the e coefficients 
in Eq. (j29p to zero); 2) ignoring the second order radiative contribution (i.e., use SL^ = and the e coefficients from 
Eqs. ([5(I1) - ([?T|) ): 3) including all second order corrections (i.e., as (2) but now with SL^ from Eq. ([5^ ). The results 
are given in Table [D From this Table we learn that the second order radiative bit which is missing from the self-force 
formalism appears to be relatively unimportant for the gravitational wave phasing, since the number of gravitational 
waveform cycles changes very little as we change the approximation used to compute them. 



-4.81022, h = 32.5926, /1.5 = -48.2687, h = -241.122; 
-4.03003, h = 18.7259, h.^ = -48.712, h = -45.5372; (32) 




M/Mq 




0.6 


1.4 


10 


Number of 


cycles with no second order corrections 


133312.1 


127501.3 


98643.4 


Number of 


cycles with no second order radiative corrections 


133311.9 


127500.8 


98641.1 


Number of 


cycles with all second order corrections 


133312.5 


127502.1 


98647.2 



TABLE I: In this table we present the number of gravitational waveform cycles that are generated during the last year of inspiral of COs 
with masses of 0.6Mq, IAMq and IOMq into a 10^ Mq Schwarzschild black hole. We consider three different approximations, as described 
in the text. 



IV. NOISE INDUCED PARAMETER ERRORS 



We now want to look at the parameter estimation accuracy which LISA observations are likely to achieve, for 
which we need waveforms in addition to the orbital phase evolution. We use the standard numerical kludge approach 
described in Babak et al. Q and refer the reader to p/7| for background. The waveform is generated by using a 
flat-space quadrupole wave generation formula applied to the trajectory of the inspiralling object in Boyer-Lindquist 
coordinates, which are identified with spherical-polar coordinates in the flat-space. In the transverse-traceless (TT) 
gauge, the metric perturbation hij{t) is given by 



h^J it) = I [p.kPjl - IPrjPk^j i"' , (33) 

where D is the distance to the source, rjij is the flat metric, the projection operator is given by = rjij — fiifij, hi 
is the unit vector in the direction of propagation, and is the second time derivative of the inertia tensor flT'l . In 
the EMRI framework it takes the form J*^' = mr^{t)r^{t), where r*(t) represents the position vector of the compact 
object with respect to the MBH in the pseudo-flat space. 

A. Implementation of EISA's response function 
Following [1^, the LISA response may be written as 
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I' 



11 t22\ „ _ rl2 



where a = /, // refers to the two independent Michelson-Hke detectors that constitute the LISA response at low 
frequencies. The functions A~^'^{t) arc the polarization coefficients given by 

A+ ^ ~a+[l + {a-nf], A"" ^ 2a-^{a ■ n), (35) 

where a is a unit vector along the SMBH's spin direction, and a+, Ox are given by 

1 

+ 2 

The antenna pattern functions F^^ are given by 

F/ = i(l + cos^ 61) cos(20)cos(2V')- cos 6* sin(20)sin(2-0), 

F]" = i(l + cos2 6l)cos(20)sin(2V') + cos6'sin(2(/))cos(2V'), (36) 



F+ = -(1 + cos^6l)sin(2(/))cos(2?/;) + cos6'cos(20)sin(2-0), 



F/^ = -(l + cos^ 61) sin(2<?!>)sin(2V^)- cos 61 cos(20)cos(2V'). (37) 

The various angles in the previous expressions represent the source's sky location in a detector based coordinate 
system, {0,(j)), and the polarization angle of the wavefront, V'- These can be re-written in a fixed, ecliptic-based 
coordinate system. If we denote the source co-latitude and azimuth angles and the direction of a in this fixed 
coordinate system by {9s,4>s) and {9k,4>k) respectively, then 

1 /3 - 

cos 9{t) = -cos9s-^sin9scos[^o + 27:it/T)-^s], 



</,(t) 




V3 cos 9s + sin 9s cos[</)o + 27r(t/T) - (j)s] 
2 sin 9s sin[^o + 2TT{t/T) - (j)s] 



V3 

tani/j = ^^cos9k — sin 9 k cos[(j)o + 2-11 (t/T) — (j)K] 

— cos 9{t) [cos 9k cos 9s + sin 9k sin 9s cos((/)a' ~ ^^'s)] j" 

{1 
- sin6'K sin 6*5 sin((/)A' ~ (j)s) ^ cos((/)o + 27rt/r) 

{cos 9k sin sin (j)s — cos sin 9k sin 0^-} 

\/3 - 1 

— sin(0o + 27rt/r) (cos05 sin6'A- cos0if — cos6'a' sin6'5 cosc/jg) >, (38) 

where 0o, ao are constant angles which represent the orbital and rotational phase of the detector at t = 0. We will set 
both of these to zero in our analysis. Additionally, T is the orbital period, which is 1 year. Barack and Cutler Q write 
these expressions in terms of 9lt4'l which specify the direction of the compact object's orbital angular momentum 
in the ecliptic-based system. In this case the orbits under consideration are circular and equatorial, so the angular 
momentum vector of the orbiting body does not precess about the SMBH's spin a and 9k = ^l, 'Pk = 4'l- 

The last ingredient in the detector response is the Doppler phase modulation. If <l>(i) denotes the phase of the 
waveform, the inclusion of the Doppler modulation shifts the phase as follows [1| 

$(t) ^ $(t) + 2^i?sin6's cos[27r(t/T) - 0s], (39) 
where R = lAU/c = 499.00478s and d^/dt is the azimuthal velocity of the orbit, cf. ([2]). 



B. Signal analysis 



In this section we will briefly review the basics of signal analysis. We can think of a GW detector as a linear system 
whose input is a GW we want to detect and whose output is a time series. This time series is a combination of both 
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the instrumental noise and a true GW signal. For LISA, the output of each of the two equivalent two arm Michelson 
detectors can be represented as 

Sa{t) = ha{t) + na{t), a = I, II. (40) 

The detection problem is to identify ha{t) in the presence of Uait). Making the usual assumption that each Fourier 
component of the noise, ha{f), is Gaussian distributed, and uncorrelated with other Fourier components (i.e., the 
noise is stationary), the ensemble average of the Fourier components of the noise have the property 

(nM) MfT) - l^if - nsM)Sc.0. (41) 

This relation defines the one-sided spectral density of the instrumental noise, Sn{f)- For LISA, the spectral density 
is the same in the two detectors. 

The problem we have to solve is how to estimate the parameters of the system. That is, if the expected waveform 
h(t; 0) depends on parameters — {9i, 9n}, we want to reconstruct the most probable value of the parameters of 
the source and compute their respective errors. To effectively address this problem we need to calculate the probability 
of the parameters given the data, i.e., the posterior probability. We may define a natural inner product on the vector 
space of signals, which for any two signals Pa{t), qa{t), takes the form 



/"OC 

(piq) ^ 2^ / mmM) +pM)cif)]/Sn{f)df. 



(42) 



The probability distribution fimction for the noise n{t) is given by 



p(no) = iVexp(-(^L (43) 



where TV is a normalization factor. Assuming we have made a detection, the output has the form s{t) — h(t; ^truc) + 
no(<), where no(i) is the specific realization of the noise and ^truo is the unknown true value of the parameters. The 
probability of measuring this output signal is then 

p (h(0true) I S) = Npo{0true)exp (^(h(0truc) I s) - ^ (h(^truc) I h(0truc)) j , (44) 

where po(^true) is a prior probability p^ . This probability distribution function encodes all the information we can 
extract from the data stream. We want to determine the most probable value of the parameters and estimate the 
likely errors in these estimates. The maximum likelihood estimators of the parameters are those that minimize the 
product (s — h I s — h). This choice, which also provides the highest SNR in a matched filtering search, is given by 

4[M^^)] = 4L^. (45) 

If the SNR is large, the parameter estimation error is likely to be small. Ignoring the prior, we can expand Eq. (|44|) 
about the peak, 6 by setting 0'^ — 0^ + and find 

ah 

Tij is the Fisher Information Matrix. For large SNR, the expectation value of the errors A^* is given by 

{Aff'A0^) = (T-^y^ + 0{SNR)-\ (47) 

In our analysis, we employ a simplified definition of the inner product, Eq. (|42p . For white noise, i.e., <S'„(/) = const., 
it would take the simple form 25*"^ J^a I^oo Pa{t)qa{t)dt, by Parseval's theorem Q. If we define the "noise-weighted" 
waveform 

we can rewrite the Fisher matrix approximately as 

Tab^2y2l daK{t)dbK{t)dt. (49) 
This approach was also used by Barack and Cutler [4j| . 



p{A0\s) = Ue-^^^^^''^<'\ T,j 



(46) 
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C. Noise model 



The function Sh{f) is the total LISA noise, which has three components: instrumental noise, confusion noise from 
short-period galactij; binaries, and confusion noise from extragalactic binaries. We use the same prescription as in 
Barack and Cutler 4] , namely 

S,.{f{t)) = mm{Sr'if)/expi-nT-Un<iN/df) + Sl-<^^\f), 

5r(/) + 5f(/) + 5r^'(/)}, (50) 
where the various components have the following analytic forms Q 

= 9.18 X 10-''^f-^ + 1.59 X lO^''^ + 9.18 x 10"^^^ Hz"\ 



-7/3 

Sfif) = 2.1 X 10-^5 ( -L- ) Hz-\ 



VlHz J 

^cx. gal ^ 4 2 X 10-4^ (^^^ ^ Hz-1. (51) 

Here dN/df is the number density of galactic white dwarf binaries per unit GW frequency, and k is the average 
number of frequency bins that are lost when each galactic binary is fitted out. We use 

diV ^ „ _i /IHz^ 



2xlO-^Hz-i( — ) , '^T^iLion - 1-5/y, (52) 



3 yr, and k = 4.5. 



D. Parameter estimation error results 



To estimate errors using the inverse Fisher Matrix we used fixed values of the intrinsic parameters of the source and 
carried out a Monte Carlo simulation over values of the extrinsic parameters. The error estimate is SNR dependent, 
but we quote results at fixed SNR= 30. We first compute the Fisher Matrix for a source at D = IGpc, and the 
corresponding SNR from the expression 




We then multiply the errors from the inverse Fisher Matrix by (SNR/30) to normalise to SNR= 30. We considered 
one year observations and chose the initial radial coordinate, pq, such that the compact object would reach the last 
stable orbit after one year of inspiral. The Fisher Matrix has ten dimensions. Four of these are intrinsic parameters, 
namely In m, InM, (7,po- The other six are extrinsic or phase parameters. We summarize the physical meaning of the 
parameters in Table [TTl 

For these Monte Carlo simulations, we evolved the orbit via Eq. (|4]), i.e., evolving rather than p, included the full 
geodesic frequency in Eq. ([6]) and included the higher order Teukolsky fits in Eq. In this way, we included lower 
order effects as accurately as possible as discussed earlier. The results of the Monte Carlo simulations are summarized 
in Tables [Till IIVI and fVl We considered "typical" systems with M — IO^A/q, q — 0.9 and three different values of 
TO = 0.6Mq, 1.4M0, IOMq to represent inspirals of white dwarfs, neutron stars and black holes, respectively. For 
comparison to the self force we also considered a case with q = 0, for which we also ignored g as a parameter in the 
Fisher Matrix. The tables list the mean, standard deviation, median and lower and upper quartiles of the distribution 
of Fisher Matrix errors computed in the Monte Carlo simulation. In Figure |4] we show sample histograms of the Fisher 
Matrix errors in the intrinsic parameters computed from the Monte Carlo simulation for the m — IOMq system. 

For the to = IOMq case we also carried out Monte Carlo simulations with (a) truncated at 2PN order, which 
did not significantly affect the results for the q = case, but made the q — 0.9 errors greater by about an order of 
magnitude; and (b) p evolved directly via the 2PN expression, Eq. ([5]), rather than evolving Lz, which made the errors 
for the q = case about an order of magnitude greater, but which did not significantly change the q = 0.9 results. 

The Tables show estimates for the noise-induced parameter errors if the last year of inspiral is observed. However, 
LISA may observe a binary system in various stages of evolution. In Figures [5] and [6] we show how the error estimates 
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In m 


mass of CO 


InM 


mass of SMBH 


<? 


magnitude of (specific) spin angular momentum of SMBH 


Po 


Initial radius of CO orbit 


4>o 


Initial phase of CO orbit 


ds 


source sky colatitude in an ecliptic-based system 


4>s 


source sky azimuth in an ecliptic-based system 


Ok 


direction of SMBH spin (colatitude) 


<t>K 


direction of SMBH spin (azimuth) 


D 


distance to source 



TABLE II: This table describes the meaning of the parameters used in our model. The angles {Os,<t>s) and {OK,<t>K) are defined in a 
fixed ecliptic-based coordinate system. 



Model 


Parameter 


log(m) 


log(A/) 


log('7) 


log(po) 


log(<?!)o) 


log(es) 


log(<;!>s) 


log(eA-) 


log (0k) 


log(D) 


Kludge 
(q=0.9) 


Mean 


-4.04 


-3.73 


-4.53 


-3.88 


-0.73 


-1.63 


-1.74 


-0.93 


-0.79 


-1.58 


St. Dev. 


0.121 


0.125 


0.129 


0.125 


0.864 


0.163 


0.205 


0.750 


0.773 


0.391 


L. Qt. 


-4.14 


-3.82 


-4.60 


-3.96 


-1.33 


-1.75 


-1.90 


-1.51 


-1.38 


-1.85 


Med. 


-4.04 


-3.72 


-4.52 


-3.87 


-1.07 


-1.62 


-1.76 


-1.16 


-1.02 


-1.72 


U. Qt. 


-3.94 


-3.63 


-4.43 


-3.77 


-0.34 


-1.50 


-1.65 


-0.53 


-0.37 


-1.41 


Self Force 
(q=0) 


Mean 


-4.37 


-4.79 


N/A 


-4.96 


-0.70 


-1.57 


-1.64 


-0.90 


-0.76 


-1.26 


St. Dev. 


0.128 


0.124 


N/A 


0.124 


0.904 


0.136 


0.204 


0.770 


0.783 


0.407 


L. Qt. 


-4.44 


-4.85 


N/A 


-5.02 


-1.33 


-1.66 


-1.78 


-1.48 


-1.35 


-1.54 


Med. 


-4.34 


-4.77 


N/A 


-4.94 


-1.05 


-1.54 


-1.66 


-1.13 


-0.95 


-1.40 


U. Qt. 


-4.27 


-4.69 


N/A 


-4.86 


-0.36 


-1.45 


-1.55 


-0.54 


-0.41 


-1.12 



TABLE HI: Summary of Monte Carlo over Fisher Matrix errors for black hole systems (m = IOMq). We show the mean, 
standard deviation, median and quartiles of the distribution of the logarithm to base ten of the error in each parameter. Results 
are given both for the kludge model (with conservative corrections to 1.5PN order) and for the self-force model, as indicated. 
The angles 4>o and ao, specifying LISA's position and orientation at t = 0, are set to zero. 



vary as a function of the time remaining until plunge at the start of the observation. These results assume a one year 
LISA observation with fixed parameters m = lOA/©, M = lO^M©, q = 0.9, $o = 0, 65 = 7r/4, (j)s ^ 0,9k = 7r/8, (j)K = 
0. We present the results in two ways, (1) the distance to the source is adjusted to keep the total SNR= 30 over 
the observation; (2) the distance to the source is fixed and chosen so that the SNR= 30 over a one year observation 
starting one year before plunge. 

Our results for black hole systems, m = IOM0, are broadly consistent with existing results in the literature [1]. 



Model 


Parameter 


log(m) 


log(M) 


log(g) 


log(po) 


log(<^o) 


log(es) 


log (0s) 


log(eK) 


log (0k) 


log(O) 


Kludge 
(q=0.9) 


Mean 


-4.60 


-3.94 


-4.59 


-4.06 


-0.75 


-1.71 


-1.88 


-0.96 


-0.82 


-0.92 


St. Dev. 


0.109 


0.118 


0.122 


0.118 


0.882 


0.199 


0.212 


0.766 


0.797 


0.399 


L. Qt. 


-4.70 


-4.02 


-4.65 


-4.14 


-1.36 


-1.87 


-2.03 


-1.56 


-1.40 


-1.19 


Med. 


-4.61 


-3.94 


-4.58 


-4.05 


-1.09 


-1.71 


-1.89 


-1.19 


-1.03 


-1.06 


U. Qt. 


-4.52 


-3.85 


-4.49 


-3.97 


-0.42 


-1.54 


-1.77 


-0.60 


-0.45 


-0.79 


Self Force 
(q=0) 


Mean 


-4.12 


-4.87 


N/A 


-5.06 


-0.76 


-1.63 


-1.71 


-0.94 


-0.80 


-0.62 


St. Dev. 


0.136 


0.130 


N/A 


0.131 


0.852 


0.155 


0.217 


0.731 


0.767 


0.379 


L. Qt. 


-4.19 


-4.95 


N/A 


-5.13 


-1.37 


-1.74 


-1.86 


-1.52 


-1.40 


-0.88 


Med. 


-4.10 


-4.86 


N/A 


-5.04 


-1.08 


-1.61 


-1.73 


-1.17 


-1.01 


-0.76 


U. Qt. 


-4.02 


-4.78 


N/A 


-4.96 


-0.41 


-1.49 


-1.62 


-0.58 


-0.40 


-0.48 



TABLE IV: As Table|llll but for neutron star inspirals (m = 1.4Mq). 
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Model 


Parameter 


log(m) 


log(Af) 


log(g) 


log(po) 


log(<^o) 


log(es) 


logics) 


log(eK) 


log(<?!>K) 


log(D) 


Kludge 
(q=0.9) 


Mean 


-4.78 


-3.93 


-4.52 


-4.03 


-0.73 


-1.74 


-1.88 


-0.94 


-0.80 


-0.57 


St. Dev. 


0.119 


0.182 


0.195 


0.182 


0.892 


0.211 


0.228 


0.781 


0.813 


0.410 


L. Qt. 


-4.82 


-4.00 


-4.57 


-4.10 


-1.36 


-1.89 


-2.04 


-1.56 


-1.43 


-0.85 


Med. 


-4.76 


-3.91 


-4.49 


-4.01 


-1.07 


-1.74 


-1.90 


-1.15 


-1.01 


-0.71 


U. Qt. 


-4.70 


-3.83 


-4.41 


-3.93 


-0.36 


-1.56 


-1.77 


-0.55 


-0.41 


-0.43 


Self Force 
(q=0) 


Mean 


-3.99 


-4.90 


N/A 


-5.08 


-0.70 


-1.64 


-1.74 


-0.90 


-0.77 


-0.26 


St. Dev. 


0.139 


0.132 


N/A 


0.133 


0.928 


0.157 


0.196 


0.799 


0.819 


0.424 


L. Qt. 


-4.05 


-4.98 


N/A 


-5.16 


-1.37 


-1.76 


-1.89 


-1.52 


-1.39 


-0.56 


Med. 


-3.96 


-4.88 


N/A 


-5.06 


-1.05 


-1.63 


-1.74 


-1.13 


-0.98 


-0.42 


U. Qt. 


-3.88 


-4.80 


N/A 


-4.98 


-0.29 


-1.51 


-1.64 


-0.48 


-0.36 


-0.11 



TABLE V: As Table |TTT1 but for white dwarf inspirals (m = O.6M0). 



We find typically that an EMRI observation can determine the CO and SMBH masses and the SMBH spin within 
fractional errors of ^ 10~^, 10~'^'^ and 10~^'^ respectively. We expect to determine the location of the somxc on 
the sky to within 10^"^ steradians, and determine the SMBH spin orientation to within lO^'^'^ steradians. However, 
these results are based on a different waveform model. Our model is based on true geodesies of the Kerr space-time, 
so among other things we have a more accurate frequency at plunge. The various other modifications, such as the 
inclusion of the conservative self-force corrections hopefully means this model is providing reliable and independent 
error estimates. One thing that is significantly different in our results is the estimate of LISA's ability to determine 
parameters for the neutron star and white dwarf systems. Our error estimates for these systems are even smaller 
than those for the black holes. The improvement arises because, by evolving rather than p, the approach to the 
innermost stable circular orbit (ISCO) shows the proper rapid change in p, which carries more information about the 
intrinsic parameters. This rapid inspiral to ISCO is missed in PN evolutions. In a sense, these results are optimistic 
as, given the event rate [3], we are unlikely to observe such systems in the last year before plunge, and the SNR 
is likely to be lower, while we have normalised all results to SNR= 30 in this paper. These results are nonetheless 
interesting. 

In the Monte Carlo simulations that used an evolution of p instead of (case (b) mentioned above), we found 
results that were close to the ones quoted here and were also in good agreement with the results obtained by Barack 
and Cutler's [1]. However, there were subtle differences. In particular, the distribution of the errors, cf., Figure SI 
is much broader when using the 2PN p evolution. In addition, the error estimate for the SMBH spin is an order of 
magnitude worse using the 2PN p evolution. 

These various results provide reassurance that our estimates from the kludge model and those by Barack & Cutler 
niay provide a fair reflection of what will be achieved in practice with LISA. Our Monte Carlo results are also 
the first of their kind for the EMRI problem to appear in the literature and illustrate the spread in errors that arise 
from randomisation of the extrinsic source parameters. The variation in the parameter determination accuracy as 
a function of time remaining until plunge indicates that the accuracy of determination of the intrinsic parameters 
depends significantly on the which part of the inspiral is observed, while the determination of the extrinsic parameters 
is equally good no matter which part of the inspiral is detected. This makes sense, as most of the orbital evolution 
takes place during the last year of inspiral and so the signature of the intrinsic parameters is most strongly present 
there. In contrast, the determination of the extrinsic parameters, such as the position in the sky and orientation of 
the SMBH spin, comes primarily from the modulations induced in the waveform by the detector motion, and so do 
not change significantly provided a full year of inspiral is observed. 

V. MODEL-INDUCED PARAMETER ERRORS 

The FM gives an estimate of the errors that arise due to noise in the detector. But that is not the only source 
of error. The kludge we have built is only approximate and hence the kludge waveform that is the best-fit to the 
data may have different parameters to the true waveform, introducing another parameter error. We shall call these 
"model" errors. Cutler and Vallisneri [9| recently introduced a framework that allows us to estimate the magnitude 
of the model errors. This relies on knowing what the "true" signal is, which we do not. However, we can use the 
formalism to estimate how important for parameter estimation the inclusion of various terms in the waveform model 
may be. In this section, we look specifically at what effect ignoring the conservative self-force corrections will have 
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FIG. 4: Distributions of errors in the intrinsic parameters as computed from the Monte Carlo simulations of the inverse Fisher 
Matrix. We show the system with m — lOM©, M = lO^M©, q — 0.9 and show the errors in ln(m), In(A/), ln{q) and In(po) 
respectively. 



on parameter determination. 

The noise error scales as 1/SNR (cf. (|47|) ). but the model errors are independent of SNR. Therefore, it might be the 
case that for the EMRI systems with highest SNR, theoretical errors could dominate the total parameter-estimation 
error. In general, if crossing a term out in the model gives rise to a parameter error that is comparable to or smaller 
than the noise-induced error, then we can safely ignore that term in our search template. We want to estimate 
theoretically what error would result from omitting conservative corrections from the waveform model, or by using a 
model with incomplete conservative corrections. So, we take the "true" waveform, /igRi to be the kludge waveform 
including all conservative pieces. We then estimate the model errors when we search for /iqr using templates, /iap, 
that include none or only part of the conservative corrections. 

We now briefly review the Cutler and Vallisneri model error formalism We consider two manifolds embedded 
in the vector spaces of data streams. One of them is covered by the true waveforms {hQ^{6^)} and the other one by 
the approximate waveforms {hAp{0^)}- Given a signal s — h.Q^{6^) + n, the best fit 6*' is determined by the condition 



Furthermore, at first order in the error A6'*(n) e 



»(n)- 



Eq. (E 



0. 



takes the form 



T-\d)) (a,hAP(0) 



hGR(^)-hAp(0) 



(54) 



(55) 



where the Fisher matrix is evaluated using the approximate waveforms ^ij{d) = {dih.Ap{0)\djhAp{d))- 

Relation ([55)1 clearly shows us that, at leading order, A9^ is the sum of two contributions. The first one is due to 
noise in the detector, A„0*, whereas the second one, Ath6'% is the contribution due to the inaccurate waveform. These 
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FIG. 5; We show plots of noise-induced errors in a one year observation as a function of time remaining until plunge at the 
start of the observation. Two cases are illustrated, namely (1) the distance to the source is normalized by means of weighting 
the noise errors by a factor of SNR/30 (solid line), (2) the distance to the source is fixed as in (1) but only for the last year 
before plunge (dot-dash line). The plots correspond to a lOM© CO orbiting around a 10® M© SMBH with parameter spin 
q = 0.9. We have set the various extrinsic parameters as follows: "l>o = 0, = 7r/4, 4>s = 0, 6k ~ tt/S, (pK ~ 0. These four 
plots show the errors in the intrinsic parameters. 



are given, respectively, by 



A„0^ = {T-'ie)] {d,hAP{e) 



n 



= r-^(^^) a,hAP(0) hGR(0)-hAP(0) 



(56) 



If we knew both 9 and the noise realization n then these equations would allow us to determine 9. However, exper- 
imentally we are only able to determine the hAp(^) that is the best fit waveform for a given data stream, s, and we 
are unsure about the error A6 = 9 — 9. In addition, we do not know 9 in Eq. (|56p . At leading order, we can replace, 
hGR(^) - hAp(^) by hGR(6') - Hap (6*), obtaining 



(r-\9)y' (a,hAp(0)|hGR(0) - hAp(0) 



(57) 



This relation is both noise and SNR independent. This property, along with the fact that Ath^* is not averaged out if 
the same event is measured by a large number of nearly identical detectors leads us to consider Ath^* as a systematic 
error. 

Cutler and Vallisneri found that the leading order approximation, Eq. ([57)1 . was not very good, unless the 
waveform is re-written in an amplitude-phase form 



(58) 
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FIG. 6: Same as in Figure [5] but for extrinsic parameters. 



The amplitude A and phase ^' are given by 

\/3 m 



Ai = 

An = 



2D \/ ' + 



A\Fl^ + AlFl,, *, = 20 + V/, 



2D 



AlFl^^ + Al Ffj^ ^ , = 20 + V// , 



(59) 
(60) 



where A+^x are given by fSS]) with A+ — Ax — 2p^(/i^, i^{/,7/;+,x} are given by ([36]) and p7)) . and the other various 
quantities are given by 



'ijjj — arctan 



AxFl^y 

'A+Fj, 



tjjii — arctan 



AxFlI^y 

A+Fjj, 



(61) 
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The first order approximation to this expression 

Ath^* « (r-i(^?))'^ ( [AA + iAA*]e^| a,hAp(0)), (62) 

at e 

was found to give rehable results when compared to more accurately computed error estimates |91] , so we use this form 
again here. Equation ([5^ behaves better than Eq. ([57]) since the difference between two waveforms, hAp(0) — hAp(^), 
is not very well approximated by the first term in its Taylor expansion. The differences in both the amplitude 
and phase of the waveform are individually well approximated by the linear terms in the Taylor series [9|. In fact, 
([57)1 is reliable only as long as the phase difference between the two waveforms is much less than one radian, i.e. 
AO^dj'^ApiO) <C 1, whereas for we just require A9^ A9^ didj"^ ap{9) ^ 1- This condition is much less restrictive 
than the former one. 

We can now use equation (|62p to estimate the magnitude of the parameter errors that arise from inaccuracies in 
the template waveform. At present, accurate waveforms including all first order conservative self-force corrections 
are not known. We want to estimate how relevant these corrections are for parameter estimation. We can use the 
kludge model for this, since we have now included conservative corrections to 2PN order and can therefore turn these 
corrections on and off. This will allow us to compute the ratio TZ of the model error that arises from omitting the 
conservative part of the self-force to the Fisher Matrix error for each of the 10 parameters in the model. This ratio 
will indicate the importance of the conservative corrections for parameter determination. If 7?, 1 then the estimates 
obtained from a model that ignores the conservative piece should still be reliable, but ii TZ >> 1 then it is clear 
that the parameter estimates would not be trustworthy. The ratio of parameter errors to FM errors obtained from a 
Monte Carlo simulation are summarized in Tables IVlllVIII and I Villi As before we have evaluated these model- induced 
parameter errors using expression ^ to evolve p, and using the full Teukolsky fit expression for L^, Eq. ([3]). We 
have also done Monte Carlo simulations with truncated at 2PN order, and evolving p at 2PN order directly using 
Eq. (O. The results were largely consistent between all three simulations. 

We quote results for the same three test systems that we used previously, and consider five different comparisons — 
for the first two we took the "true" waveform to be our kludge waveform with 2PN conservative corrections, and took 
the template to be a kludge waveform with either no conservative corrections ( "OPN" ) or with conservative corrections 
to 1.5PN order ("1.5PN"); for the latter three comparisons we used the self-force model in the three varieties listed 
in Table H]— (1) SL^ = and all = in Eq. ^ ("1st order"); (2) SL^ = and the e coefficients set to the 
values in Eqs. ([30)) - ([3T|) ("incomplete"); (3) as (2) but now with 6Lz from Eq. ([32]) ("2nd order") — and did pairwise 
comparisons. 

In Figure [71 we also show the full distribution of the error ratio, TZ, for the intrinsic parameters, as computed from 
the Monte Carlo simulation. These are shown for the black hole systems in the kludge 2PN vs 1.5PN comparison. 
These plots show that, in this particular case, the vast majority of sources fulfill the condition that the model errors 
are smaller than the noise induced errors. For those binary systems that do not fulfill this condition, the upper bound 
on the error ratio 7^ is 7^ < 3. Tables IVll IVlIl and IVlIIl indicate that, for both kludge model comparisons, the ratio 
TZ < 1 for all parameters at most points in the parameter space for all three types of inspiral. In fact, for any of the 
model parameters, less than 0.15% of points in the Monte Carlo runs had 2 < TZ < 3. These results suggest that 
including conservative corrections is not essential for accurate parameter determination, but including them up to 
2PN order will certainly reduce the model errors further. 

Things are not so clear cut for the self-force comparisons, however. For the neutron star and white dwarf inspirals, 
the ratio TZ < 1 for practically every element of the parameter space. In fact, for any parameter, less than 0.1% of 
points of the Monte Carlo runs satisfy 2 < 7?, < 3. For black holes, the model errors appear generally larger, with 
ratios typically greater than 1 in the comparisons to the "2nd order" model that was built from the kludge. One 
thing that must be born in mind for these comparisons is that the kludge model was constructed by comparison 
to a weak-field post-Newtonian expansion, but we are now comparing it to fully accurate self-force computations in 
the strong-field. Thus, the reason for the apparently larger discrepancy may be that the kludge model is itself not 
accurate enough in the regime where we are comparing. This was the reason that we also did the "incomplete" vs 
"1st order" comparisons for the self-force models. These expressions do not use the kludge model, but are in the spirit 
of the other kludge comparisons in that we are crossing out the last term in the model and looking at what effect this 
has. Those results show generally smaller errors, with more than 50% of points having TZ < 1. The conclusion is that, 
in the worst case, a template that omits conservative corrections will identify parameters that are (conservatively) 
within 10 Fisher Matrix errors of the true parameters, and this region of parameter space can then be followed up 
using more accurate waveforms if these are available. We also conclude that the kludge model with 2PN conservative 
corrections is a reasonably good approximation to the self-force model, since although the ratios for those comparisons 
are larger, they are still manageably small. 
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Model 


Parameter 


log(m) 


log(M) 


log(g) 


log(po) 


log(<?!'o) 




log(<^s) 


log(eA') 


log(<?f>K) 


log(D) 


Kludge 
2PN vs OPN 


Mean 


-0.22 


-0.23 


-0.17 


-0.22 


-0.09 


-0.20 


-0.20 


0.00 


0.03 


0.60 


St. Dev. 


0.584 


0.626 


0.572 


0.611 


0.736 


0.581 


0.534 


0.750 


0.745 


0.832 


L. Qt. 


-0.59 


-0.53 


-0.48 


-0.54 


-0.41 


-0.43 


-0.47 


-0.38 


-0.29 


0.19 


Med. 


-0.11 


-0.10 


-0.06 


-0.10 


0.01 


-0.10 


-0.09 


0.11 


0.16 


0.75 


U. Qt. 


0.25 


0.23 


0.25 


0.23 


0.39 


0.20 


0.17 


0.55 


0.54 


1.20 


Kludge 
2PN vs 1.5PN 


Mean 


-0.45 


-0.43 


-0.38 


-0.43 


-0.33 


-0.40 


-0.44 


-0.23 


-0.20 


0.36 


St. Dev. 


0.598 


0.591 


0.540 


0.594 


0.744 


0.578 


0.606 


0.720 


0.733 


0.841 


L. Qt. 


-0.75 


-0.72 


-0.65 


-0.71 


-0.58 


-0.63 


-0.66 


-0.53 


-0.53 


0.03 


Med. 


-0.32 


-0.31 


-0.29 


-0.32 


-0.20 


-0.28 


-0.30 


-0.09 


-0.08 


0.54 


U. Qt. 


-0.01 


-0.00 


0.01 


-0.00 


0.12 


-0.01 


-0.01 


0.24 


0.29 


0.91 


Self Force 
"incomplete" vs 
"1st order" 


Mean 


-0.01 


0.06 


N/A 


0.05 


-0.04 


-0.14 


-0.18 


0.11 


0.15 


0.71 


St. Dev. 


0.602 


0.571 


N/A 


0.577 


0.940 


0.814 


0.845 


0.968 


0.947 


1.010 


L. Qt. 


-0.39 


-0.31 


N/A 


-0.31 


-0.40 


-0.42 


-0.51 


-0.29 


-0.21 


0.27 


Med. 


0.07 


0.09 


N/A 


0.09 


0.07 


0.05 


0.05 


0.32 


0.32 


0.88 


U. Qt. 


0.46 


0.52 


N/A 


0.51 


0.59 


0.38 


0.41 


0.76 


0.77 


1.38 


Self Force 
"2nd order" vs 
"1st order" 


Mean 


0.12 


-0.03 


N/A 


-0.04 


-0.04 


-0.03 


-0.02 


0.14 


0.16 


0.71 


St. Dev. 


0.456 


0.651 


N/A 


0.663 


0.740 


0.510 


0.522 


0.734 


0.711 


0.784 


L. Qt. 


-0.15 


-0.45 


N/A 


-0.46 


-0.40 


-0.34 


-0.33 


-0.26 


-0.25 


0.25 


Med. 


0.16 


0.08 


N/A 


0.08 


0.01 


0.01 


0.10 


0.21 


0.24 


0.81 


U. Qt. 


0.43 


0.50 


N/A 


0.49 


0.44 


0.34 


0.38 


0.67 


0.67 


1.31 


Self Force 
"2nd order" vs 
"incomplete" 


Mean 


0.34 


0.28 


N/A 


0.27 


0.29 


0.26 


0.22 


0.44 


0.45 


1.01 


St. Dev. 


0.538 


0.630 


N/A 


0.632 


0.794 


0.555 


0.605 


0.758 


0.770 


0.857 


L. Qt. 


-0.00 


-0.17 


N/A 


-0.16 


-0.14 


-0.06 


-0.12 


-0.03 


0.05 


0.53 


Med. 


0.38 


0.35 


N/A 


0.35 


0.36 


0.35 


0.36 


0.54 


0.56 


1.15 


U. Qt. 


0.74 


0.77 


N/A 


0.76 


0.85 


0.68 


0.67 


0.99 


1.01 


1.64 



TABLE VI: Summary of Monte Carlo simulation results for the ratio of model errors to Fisher Matrix errors for black hole 
systems (m — IOAIq) using the Teukolsky fit expression for dL^/dt. We show the mean, standard deviation, median and 
quartiles of the distribution of the logarithm to base ten of the ratio for each parameter. Results are given for various 
comparisons, as indicated and described in the text. A comparison "A vs B" uses model A as the true waveform, and model 
B as the search template. 

The last two rows of the tables also indicate that it may be better to work consistently to a certain order in mass 
ratio, rather than to include an incomplete term at higher order. The "incomplete" model includes second order in rj 
corrections to the evolution of the orbital radius that arise from the first order conservative part of the self-force, but 
not those that arise from the second order radiative part. The "1st order" model includes no second order corrections 
to the rate of change of radius. The tables indicate that the "1st order" model actually leads to smaller parameter 
errors (when searching for the kludge waveform) than the "incomplete" model, although the difference between the 
two cases is relatively small. 

One other thing that we should point out is that wc are normalising to SNR=: 30 for all the results. This is an 
estimate of the threshold SNR that will be required for EMRI detection but nearby events may have SNR as 
high as several hundred [13] ■ For such systems, the noise induced errors will be smaller than those quoted here by a 
factor of SNR/30, while the model induced errors will be the same. Thus, for accurate parameter determination for 
the loudest sources LISA detects, these results suggest that it will almost certainly be necessary to include first-order 
conservative and, if possible, second-order radiative terms in the model. 

In Figures [8] and [9] we show how the model and FM errors vary for a one year observation as a function of the 
time remaining until plunge at the start of the observation. This is for a black hole inspiral, m = lOM©, with the 
following fixed values of the extrinsic parameters, $0 = 0,9s = '^/4:,(t>s — 0,0k — ''^l^,<t>K — 0. In each plot, the 
noise errors are weighted by a factor SNR/30 to ensure the SNR is 30 over the one year observation. We have also 
produced similar plots with the source kept at a fixed distance, but we do not show these here as they do not look 
very different. 



Model 


Parameter 


log(m) 


log(M) 


log(g) 


log(po) 


log(<?!'o) 




log(<^s) 


log(eA') 


log(<?f>K) 


log(D) 


Kludge 
2PN vs OPN 


Mean 


-0.43 


-0.35 


-0.33 


-0.36 


-0.35 


-0.33 


-0.36 


-0.23 


-0.21 


0.37 


St. Dev. 


0.563 


0.587 


0.586 


0.626 


0.728 


0.589 


0.601 


0.728 


0.754 


0.832 


L. Qt. 


-0.75 


-0.62 


-0.60 


-0.62 


-0.59 


-0.55 


-0.58 


-0.48 


-0.44 


0.05 


Med. 


-0.30 


-0.20 


-0.19 


-0.20 


-0.19 


-0.18 


-0.20 


-0.10 


-0.04 


0.57 


U. Qt. 


-0.00 


0.06 


0.08 


0.06 


0.08 


0.07 


0.03 


0.23 


0.23 


0.89 


Kludge 
2PN vs 1.5PN 


Mean 


-0.80 


-0.67 


-0.68 


-0.68 


-0.73 


-0.68 


-0.71 


-0.57 


-0.56 


0.02 


St. Dev. 


0.576 


0.551 


0.657 


0.573 


0.744 


0.598 


0.623 


0.711 


0.760 


0.817 


L. Qt. 


-1.09 


-0.90 


-0.88 


-0.90 


-0.98 


-0.90 


-0.91 


-0.83 


-0.77 


-0.23 


Med. 


-0.66 


-0.53 


-0.50 


-0.54 


-0.57 


-0.51 


-0.56 


-0.43 


-0.40 


0.21 


U. Qt. 


-0.37 


-0.29 


-0.28 


-0.29 


-0.31 


-0.26 


-0.33 


-0.15 


-0.12 


0.52 


Self Force 
"incomplete" vs 
"1st order" 


Mean 


-0.85 


-0.80 


N/A 


-0.81 


-0.91 


-0.96 


-0.94 


-0.77 


-0.74 


-0.17 


St. Dev. 


0.569 


0.590 


N/A 


0.584 


0.865 


0.686 


0.723 


0.870 


0.886 


0.916 


L. Qt. 


-1.18 


-1.14 


N/A 


-1.13 


-1.21 


-1.18 


-1.19 


-1.06 


-1.04 


-0.48 


Med. 


-0.77 


-0.75 


N/A 


-0.77 


-0.79 


-0.80 


-0.75 


-0.62 


-0.60 


-0.01 


U. Qt. 


-0.43 


-0.35 


N/A 


-0.36 


-0.42 


-0.51 


-0.48 


-0.27 


-0.24 


0.39 


Self Force 
"2nd order" vs 
"1st order" 


Mean 


-0.96 


-1.00 


N/A 


-1.01 


-1.24 


-1.25 


-1.21 


-1.03 


-1.03 


-0.46 


St. Dev. 


0.398 


0.420 


N/A 


0.433 


0.734 


0.576 


0.586 


0.715 


0.725 


0.851 


L. Qt. 


-1.10 


-1.14 


N/A 


-1.14 


-1.50 


-1.49 


-1.47 


-1.28 


-1.32 


-0.78 


Med. 


-0.87 


-0.92 


N/A 


-0.93 


-1.09 


-1.08 


-1.07 


-0.91 


-0.89 


-0.26 


U. Qt. 


-0.70 


-0.74 


N/A 


-0.74 


-0.79 


-0.86 


-0.79 


-0.58 


-0.57 


0.09 


Self Force 
"2nd order" vs 
"incomplete" 


Mean 


-0.60 


-0.71 


N/A 


-0.72 


-0.74 


-0.74 


-0.74 


-0.57 


-0.55 


0.03 


St. Dev. 


0.461 


0.610 


N/A 


0.611 


0.721 


0.568 


0.612 


0.747 


0.757 


0.824 


L. Qt. 


-0.85 


-1.02 


N/A 


-1.02 


-1.06 


-0.93 


-1.02 


-0.87 


-0.85 


-0.31 


Med. 


-0.54 


-0.64 


N/A 


-0.65 


-0.64 


-0.63 


-0.61 


-0.43 


-0.42 


0.16 


U. Qt. 


-0.30 


-0.24 


N/A 


-0.25 


-0.28 


-0.35 


-0.30 


-0.11 


-0.06 


0.55 



TABLE VII: As Tablc|Vll but for neutron star inspirals (m = I.4M0). 



Model 


Parameter 


log(m) 


log(M) 


log(g) 


log(po) 


log(<?!'o) 


log(es) 


log(<^s) 


log(eK) 


log(<?f>K) 


log(D) 


Kludge 
2PN vs OPN 


Mean 


-0.37 


-0.52 


-0.52 


-0.53 


-0.55 


-0.49 


-0.52 


-0.43 


-0.39 


0.17 


St. Dev. 


0.526 


0.589 


0.661 


0.614 


0.782 


0.657 


0.643 


0.786 


0.772 


0.878 


L. Qt. 


-0.57 


-0.72 


-0.70 


-0.73 


-0.78 


-0.69 


-0.71 


-0.65 


-0.60 


-0.07 


Med. 


-0.22 


-0.38 


-0.36 


-0.38 


-0.38 


-0.31 


-0.35 


-0.25 


-0.22 


0.41 


U. Qt. 


0.00 


-0.13 


-0.12 


-0.13 


-0.11 


-0.07 


-0.12 


0.03 


0.06 


0.71 


Kludge 
2PN vs 1.5PN 


Mean 


-0.71 


-0.85 


-0.86 


-0.86 


-0.96 


-0.87 


-0.86 


-0.81 


-0.79 


-0.22 


St. Dev. 


0.491 


0.518 


0.589 


0.530 


0.781 


0.632 


0.628 


0.790 


0.769 


0.870 


L. Qt. 


-0.88 


-1.06 


-1.05 


-1.06 


-1.20 


-1.09 


-1.00 


-1.04 


-1.00 


-0.48 


Med. 


-0.59 


-0.76 


-0.74 


-0.76 


-0.79 


-0.70 


-0.71 


-0.65 


-0.63 


-0.01 


U. Qt. 


-0.38 


-0.49 


-0.48 


-0.49 


-0.51 


-0.46 


-0.48 


-0.33 


-0.34 


0.32 


Self Force 
"incomplete" vs 
"1st order" 


Mean 


-1.20 


-1.17 


N/A 


-1.18 


-1.28 


-1.29 


-1.24 


-1.16 


-1.11 


-0.53 


St. Dev. 


0.540 


0.575 


N/A 


0.590 


0.757 


0.601 


0.622 


0.807 


0.728 


0.791 


L. Qt. 


-1.48 


-1.54 


N/A 


-1.54 


-1.58 


-1.49 


-1.49 


-1.42 


-1.42 


-0.87 


Med. 


-1.15 


-1.10 


N/A 


-1.12 


-1.17 


-1.16 


-1.08 


-1.05 


-1.00 


-0.41 


U. Qt. 


-0.80 


-0.70 


N/A 


-0.71 


-0.84 


-0.90 


-0.82 


-0.69 


-0.66 


-0.02 


Self Force 
"2nd order" vs 
"1st order" 


Mean 


-1.66 


-1.65 


N/A 


-1.64 


-2.15 


-2.08 


-2.09 


-2.06 


-2.01 


-1.46 


St. Dev. 


0.479 


0.472 


N/A 


0.477 


0.590 


0.518 


0.559 


0.618 


0.597 


0.730 


L. Qt. 


-1.88 


-1.88 


N/A 


-1.87 


-2.39 


-2.29 


-2.34 


-2.33 


-2.25 


-1.76 


Med. 


-1.49 


-1.47 


N/A 


-1.46 


-2.00 


-1.96 


-1.97 


-1.93 


-1.88 


-1.29 


U. Qt. 


-1.34 


-1.32 


N/A 


-1.32 


-1.77 


-1.71 


-1.71 


-1.64 


-1.62 


-0.95 


Self Force 
"2nd order" vs 
"incomplete" 


Mean 


-1.09 


-1.17 


N/A 


-1.18 


-1.24 


-1.22 


-1.19 


-1.06 


-1.01 


-0.46 


St. Dev. 


0.466 


0.581 


N/A 


0.594 


0.760 


0.572 


0.611 


0.733 


0.722 


0.834 


L. Qt. 


-1.33 


-1.48 


N/A 


-1.47 


-1.51 


-1.42 


-1.43 


-1.35 


-1.29 


-0.79 


Med. 


-1.01 


-1.09 


N/A 


-1.09 


-1.10 


-1.09 


-1.03 


-0.92 


-0.87 


-0.32 


U. Qt. 


-0.78 


-0.74 


N/A 


-0.76 


-0.79 


-0.83 


-0.78 


-0.59 


-0.55 


0.07 



TABLE VIII: As Table |Vll but for white dwarf inspirals (m = O.6M0). 
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Ratio of model error to FM error for CO's mass 



Ratio of model error to FM error for SMBH mass 



0.03 - 
0.025 - 
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Ratio of mociei error to FM error for SMBH spin 
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0.045 - 
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0.005 - 
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-3.5 -3 -2.5 -2 -1.5 -1 -0.5 

Ratio of model error to FM error 
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FIG. 7: Distribution of the ratio of the model induced error to the noise induced FM error for the intrinsic parameters, 
data is for the black hole inspiral system (m = WMq) for the kludge 2PN vs 1.5PN comparison. 



This 
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years years 

FIG. 8: We show how the noise induced and model errors vary in a one year observation as a function of the time remaining 

until plunge at the start of the observation. We show noise induced errors (solid lines), and model errors for two cases, namely 
(/iGR,/!-Ap)= (2PN, OPN) (dashed lines), and {hGR,hAp)= (2PN, 1.5PN) (dot-dash lines). This is for a lOM© CO inspiralling 
into a 10® M0 SMBH with parameter spin q = 0.9. We have set the various extrinsic parameters as follows: $0 = 0, = 7r/4, 
<ps = 0, 9k = 7r/8, 4>K = 0. The noise induced errors have been re-normalised to a signal to noise ratio of 30 over the year of 
observation. These four plots show the errors in the intrinsic parameters of the system. 
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FIG. 9: As Figure |8] but now we show the errors in the extrinsic parameters. 
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From Figures [8] and [9] we see that, for this typical source, (1) the noise induced errors in the intrinsic parameters 
increase quickly if we observe only an earlier part of the inspiral, while the model errors only increase by a small 
amount. Hence, we would only expect to be model-error dominated when observing the last year if at all; (2) noise and 
model errors for the extrinsic parameters do not change substantially as we vary the initial time; (3) the theoretical 
errors decrease as we include more conservative corrections. 

It is worth pointing out that the results quoted in Tables IVIi IVII) and IVIIII were obtained from Fisher Matrices 
which were nicely convergent over several orders of magnitude in the offset used to compute the numerical waveform 
derivatives. The noise-induced and model errors obtained from the inverses of these matrices also exhibited good 
convergent behaviour. This is a relevant observation as the matrices we are dealing with have very large condition 
numbers — the ratio of the largest to smallest eigenvalues. Typically, we found that using an LU decomposition, 
i.e., writing the FM as the product of a lower triangular matrix and an upper triangular matrix, the inverse Fisher 
matrices converged to <^ 8% over three orders of magnitude in the numerical offsets. 

VI. CONCLUSIONS 

We have constructed an improved kludge model of gravitational wave emission from circular-equatorial EMRIs 
by including conservative self- force corrections up to 2PN order. We have obtained a 2PN expression for (l which 
includes both the conservative self-force at first order in ry, and quadratic terms of the spin parameter q at 2PN order. 
Previous expressions included either conservative corrections (lH or quadratic terms of the spin parameter (l6j . but 
not both. Our model has shown us that the inclusion of conservative corrections has a relatively small impact on the 
waveform phasing, and so are not essential for source detection, but it will be useful to include them for parameter 
estimation. At this PN order, the model provides parameter determination accuracy estimates for black hole inspiral 
systems, m = IOMq, that are broadly consistent with previous results by Barack and Cutler For a typical source 
at SNR of 30, we find that a LISA EMRI observation should be able to determine the CO and SMBH masses and 
the SMBH spin magnitude to within fractional errors of ~ 10~^, 10~'^'^ and lO"'''^ respectively. We should also be 
able to determine the location of the source on the sky to within 10"'^ steradians, and determine the SMBH spin 
orientation to within lO^'''^ steradians. Our model should be more reliable in the strong field regime as we have built 
it using a true Kerr geodesic, have included conservative corrections in a physically consistent way and have evolved 
the geodesic parameters using the best available radiative flux, Lz- The fact that we get such consistent results 
provides reassurance that both estimates will be a reasonable reflection of the precision that LISA will ultimately 
achieve. We have also obtained new results for neutron star and white dwarf inspirals that indicate LISA will also be 
able to return highly accurate parameters for these systems, provided that the last year of inspiral is observed. 

We have also studied in detail the importance of the first order conservative part and the second-order radiative 
part of the self-force for parameter estimation accuracy using the formalism of Cutler and Vallisneri We have 
found that, for these sources, the model errors that arise from omitting these self-force terms are generally smaller 
than the parameter errors that arise from instrumental noise when the source has SNR= 30. In our Monte Carlo 
simulations, no points had model-error to noise-error ratios, TZ, greater than 3 and less than 0.15% of the points in 
the Monte Carlo runs lay in the range 2 < 7?. < 3. 

We have also compared our results to recently published self-force calculations, that include all first order terms, 
but nothing at higher order. Comparing to these results allowed us to assess the relative importance of the first order 
conservative and second order radiative parts of the self-force, as these affect the orbital evolution in the same way. 
These clearly showed that the missing terms were not necessary for accurate parameter estimation for the inspirals 
of neutron stars or white dwarfs. However, the results for the inspirals of black holes were less conclusive — model 
errors were typically a few times the expected parameter errors from instrumental noise. Part of the reason for this 
may have been that we derived the kludge corrections by comparison to post-Newtonian expressions in the weak field 
and then used them to test the self-force results in the strong field. Comparing the self-force waveform to a truncated 
version of the same did indicate slightly smaller error ratios, although these were still somewhat bigger than in the 
calculations based entirely on the kludge waveforms. 

Our results are the first attempt to assess the necessity of including conservative corrections in templates for 
parameter estimation with LISA. The results are not absolutely conclusive since they show neither that the model 
errors are always completely negligible nor that the model errors always overwhelm the errors from instrumental 
noise. Instead, we find that the two errors are generally comparable. This suggests that search templates can 
certainly ignore conservative corrections, but it may be necessary to follow-up with more accurate templates to get 
more precise parameter estimates. The potential problem is that our results suggest the second order radiative part 
of the self-force may be as important as the first order conservative piece and it is not clear whether templates 
that include both will be available on the timescale of LISA. Nonetheless, for all the cases we have considered, the 
model errors are only a few times the noise errors which means that a) we will obtain good estimates of the source 
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parameters, although quoted error bounds must allow for the model error; b) if more accurate templates are available, 
we will only need to use these (presumably computationally more expensive templates) in an area of parameter space 
approximately ten times the error box predicted by the Fisher Matrix. 

These results are not the full story, since we have considered a rather special class of EMRIs, namely those that 
are both circular and equatorial. If non-standard EMRI channels operate efficiently, these may make up a significant 
fraction of LISA events However, we expect orbits generically to be both eccentric and inclined to the equatorial 
plane. Although we hope our results are representative of the general case, this needs to be properly explored in the 
future. It is relatively straightforward to adapt the formalism presented here to the generic case. In general, rather 
than having just one frequency and frequency derivative, an orbit can be characterised by three frequencies and their 
derivatives (the perihelion and orbital plane precession frequencies in addition to the orbital frequency) and described 
by three parameters, e.g., a semi-latus rectum, eccentricity and inclination. Comparison of these six quantities between 
the kludge and post-Newtonian formalisms will allow determination of the coordinate transformation that relates the 
two systems and also the computation of the conservative corrections to the three frequencies. This is work that 
should be done in the near future, to verify the conclusions of this paper in a more general context. 
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